##########################################
## "Voters get what they want"          ##
##########################################
## Heinrich, Kobayashi, Long            ##
##########################################

##          
## Analysis of January survey experiment
## Script 12
## June 29, 2017


## Load and prep data from survey experiment (January)
######################################################
data <- read.csv2("data/MTurk data from February 2015.csv", quote="",sep="|")[, c("p4_screener_failures",
                                                                                  "p8_support", "p8_cost", "p8_nat_sec_just",
                                                                                  "p8_time")]
data$Time <- as.numeric(as.character(data$p8_time))
data$Failures <- data$p4_screener_failures
data$Rating <- data$p8_support
data$Just <- 0
data$Just[grepl(x=data$p8_nat_sec_just, pattern="torturing")] <- 1
data$Just[grepl(x=data$p8_nat_sec_just, pattern="counter")] <- 2
data$Costs <- data$p8_cost
data <- na.omit(data)
data <- data[, c("Rating", "Costs", "Just", "Time", "Failures")]

## Justification and Costs treatments are balanced
table(data[, c("Just", "Costs")])


## Remove people with excessive screener failures/ times
########################################################
## 1,738 observations; 1,676 retained
n <- nrow(data)
data <- subset(data, Time > 5 & Time < 180)
data <- subset(data, Failures <= 2)
n_dropped <- n - nrow(data)


## Likely broken randomization
##############################
table(data$Rating) 
0.5 * nrow(data) * 0.25
sum(data$Rating == 4) - 0.5 * nrow(data) * 0.25
## Recode
data$Rating[data$Rating == 4] <- 3


## Nonparametric bootstrap estimation
#####################################
nbs <- 20000
## Rearrange data for easier looping/bootstrapping
dat_ns <- list(subset(data, Just == 0),
               subset(data, Just == 1),
               subset(data, Just == 2))

out <- array(NA, dim=c(nbs, 3, 3), dimnames = list(Draw=NULL, Just=c("No justification", "National security",
                                                                   "National security/ CT"), 
                                                   Outcome=c("Strongly oppose", "Mildly oppose", "Any Support")))
for(i in 1:nbs)
{
  for(k in 1:3)
  {
    ind <- sample(1:nrow(dat_ns[[k]]), size=nrow(dat_ns[[k]]), replace=TRUE)
    pr1 <- 1/0.5 * (mean(dat_ns[[k]]$Rating[ind] == 1) - (1 - 0.5) * 0.25)
    pr2 <- 1/0.5 * (mean(dat_ns[[k]]$Rating[ind] == 2) - (1 - 0.5) * 0.25)
    pr3 <- 1/0.5 * (mean(dat_ns[[k]]$Rating[ind] == 3) - (1 - 0.5) * 0.50)
    out[i, k, ] <- c(pr1, pr2, pr3)
  }
}

results <- adply(.data=out, .margins=c(2,3), .fun=function(x) data.frame(q_m=median(x), q_l=quantile(x, .025),
                                                                         q_u=quantile(x, .975)))
results$Outcome <- factor(results$Outcome, levels=c("Strongly oppose", "Mildly oppose",
                                                    "Any Support"))
results$Just <- factor(results$Just, levels=c("No justification", "National security", 
                                              "National security/ CT"))


## Make graph
g <- ggplot(data=results, aes(x=Outcome, y=q_m, ymin=q_l, ymax=q_u,
                          group=Just, colour=Just))
g <- g + geom_hline(yintercept=0, size=.4)
g <- g + geom_linerange(position=position_dodge(width=.25), size=.9)
g <- g + geom_point(position=position_dodge(width=.25), size=3)
g <- g + xlab("") + theme_bw() 
g <- g + scale_colour_manual(values=c("black", "grey50", "grey80"))
g <- g + ylab("Probability") 
g <- g + guides(colour=guide_legend(title="Justification\nTreatment"))
g <- g + theme(axis.text = element_text(size=rel(1)),
               panel.grid.major.y = element_line(size=0.25, colour="grey80", 
                                                 linetype="dashed"),
               panel.grid.major.x = element_blank(),               
               panel.grid.minor = element_blank(),
               axis.ticks = element_blank(),
               strip.text = element_text(size=rel(1.2), hjust=0),
               strip.background = element_blank(),
               plot.title = element_text(size=rel(1.45), face="bold", hjust=0))
g <- g + ggtitle("Support for aid to human rights violator (January survey)")
ggsave(file="output/figures/A2-January.pdf", width=11, height=5)






